/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Patrick van der Velde
 * Marcus Cuda
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
using System;
using dnAnalytics.Math;
using dnAnalytics.Properties;

namespace dnAnalytics.LinearAlgebra.Decomposition
{
    internal abstract class AbstractLU
    {
        protected readonly Matrix mMatrix;
        private bool mComputed;
        private double mDeterminant = Double.MinValue;
        private bool mIsSingular;
        protected int[] mPivots;

        protected AbstractLU(Matrix matrix)
        {
            mMatrix = matrix.Clone();
        }

        public bool IsSingular()
        {
            Compute();
            return mIsSingular;
        }

        public Matrix LowerFactor()
        {
            Compute();
            Matrix result = mMatrix.GetLowerTriangle();
            for (int i = 0; i < result.Rows; i++)
            {
                result.ValueAt(i, i, 1);
            }
            return result;
        }

        public Matrix UpperFactor()
        {
            Compute();
            return mMatrix.GetUpperTriangle();
        }

        public double Determinant()
        {
            Compute();
            if (mIsSingular)
            {
                return 0;
            }

            if (mDeterminant == double.MinValue)
            {
                lock (mMatrix)
                {
                    mDeterminant = 1.0;
                    for (int j = 0; j < mMatrix.Rows; j++)
                    {
                        if (mPivots[j] != j)
                        {
                            mDeterminant = -mDeterminant*mMatrix.ValueAt(j, j);
                        }
                        else
                        {
                            mDeterminant *= mMatrix.ValueAt(j, j);
                        }
                    }
                }
            }
            return mDeterminant;
        }

        public int[] Pivots()
        {
            Compute();
            int[] ret = new int[mPivots.Length];
            Buffer.BlockCopy(mPivots, 0, ret, 0, mPivots.Length*Constants.SizeOfInt);
            return ret;
        }

        public Matrix Solve(Matrix input)
        {
            Matrix result = mMatrix.CreateMatrix(mMatrix.Columns, input.Columns);
            Solve(input, result);
            return result;
        }

        public void Solve(Matrix input, Matrix result)
        {
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (mMatrix.Rows != input.Rows)
            {
                throw new NotConformableException("input", Resources.ParameterNotConformable);
            }

            if (mMatrix.Columns != result.Rows)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (input.Columns != result.Columns)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (IsSingular())
            {
                throw new SingularMatrixException();
            }

            input.CopyTo(result);
            Solve(mMatrix, mPivots, result);
        }

        public Vector Solve(Vector input)
        {
            Vector result = mMatrix.CreateVector(mMatrix.Columns);
            Solve(input, result);
            return result;
        }

        public void Solve(Vector input, Vector result)
        {
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }
            if (result == null)
            {
                throw new ArgumentNullException("result");
            }

            if (mMatrix.Rows != input.Count)
            {
                throw new NotConformableException("input", Resources.ParameterNotConformable);
            }

            if (mMatrix.Columns != result.Count)
            {
                throw new NotConformableException("result", Resources.ParameterNotConformable);
            }

            if (IsSingular())
            {
                throw new SingularMatrixException();
            }

            input.CopyTo(result);
            Solve(mMatrix, mPivots, result);
        }

        public Matrix Inverse()
        {
            if (IsSingular())
            {
                throw new SingularMatrixException();
            }

            Compute();
            return ComputeInverse(mMatrix, mPivots);
        }

        protected abstract void DoCompute(Matrix matrix, int[] pivots);

        protected abstract Matrix ComputeInverse(Matrix matrix, int[] pivots);

        protected abstract void Solve(Matrix factor, int[] pivots, Matrix result);

        protected abstract void Solve(Matrix factor, int[] pivots, Vector result);

        private void Compute()
        {
            if (!mComputed)
            {
                mPivots = new int[mMatrix.Rows];

                DoCompute(mMatrix, mPivots);
                for (int j = 0; j < mMatrix.Rows; j++)
                {
                    if (mMatrix.ValueAt(j, j) == 0.0)
                    {
                        mIsSingular = true;
                        break;
                    }
                }
                mComputed = true;
            }
        }
    }
}
